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Type la supernovae result when carbon-oxygen white dwarfs in binary systems 
accrete mass from companion stars, reach a critical mass, and explode. The near 
uniformity of their light curves makes these supernova good standard candles for 
measuring cosmic expansion 1 ' 2 ' 3 ' 4 , but a correction must be applied to account for 
the fact that the brighter supernovae have broader light curves 5 . One-dimensional 
modelling, with a certain choice of parameters, can reproduce this general trend in 
the width-luminosity relation 6 ' 7 ' 8 , but the processes of ignition and detonation have 
recently been shown to be intrinsically asymmetric 9 ' 10 ' 11 ' 12 ' 13 . Here we report on 
multi-dimensional modelling of the explosion physics and radiative transfer that 
reveals that the breaking of spherical symmetry is a critical factor in determining 
both the width luminosity relation and the observed scatter about it. The deviation 
from sphericity can also explain the finite polarization detected in the light from 
some supernovae 14 . The slope and normalization of the width-luminosity relation 
has a weak dependence on certain properties of the white dwarf progenitor, in 
particular the trace abundances of elements other than carbon and oxygen. Failing 
to correct for this effect could lead to systematic overestimates of up to 2% in the 
distance to remote supernovae. 



In the most established model for SNe la, a carbon fusion flame is ignited near the 
center of the white dwarf and initially burns in a sub-sonic and turbulent deflagration, 
then transitions near the white dwarf surface to a supersonic detonation 15 ' 16 . The 
detonation is needed in order to match the observed energetics and nucleosynthesis 17 . 
The energetic stage of the explosion lasts only a second or so, but synthesizes radioactive 
56 Ni that powers the subsequent light curve and determines its luminosity. In previous 1- 
dimensional models, the 56 Ni yield depended on the choices of parameters representing 
the speed of the subsonic burning front and the density at which the front makes a 
transition to a detonation wave . However, these parameters are actually a consequence 
of multi-dimensional instabilities not captured in ID, and so were not highly constrained 
by physics. Ignition in the ID models also occurred at the very center of the star, and the 
transition to detonation happened simultaneously on a symmetric spherical shell. Both of 
these assumptions are now doubted. 

The starting point of our simulations is a standard 1.38 solar mass white dwarf composed 
of 50% carbon and 50% oxygen with a central density of 2.9 x 10 g cm" . The 
propagation of the burning is followed in 2 spatial dimensions using a level set to track 
the flame's location and a turbulent subgrid model adopted from the chemical 
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combustion community to describe its motion . The resulting debris are then post- 
processed using a multi-dimensional radiative transport code 19 to determine the emergent 
radiation. Model variations consist of how the white dwarf is ignited and the criteria for 
making the transition to detonation, both of which are based on insights from recent 3-D 
studies, as described below. Variations in the trace abundance of elements other than 
carbon/oxygen (the metallicity) of the progenitor star over a range from 1/3 to 3 times 
solar are also explored to mimic the evolution in supernova environment that may have 
occurred over aeons of cosmic time. 



Simulations show 9 ' 20 that just prior to ignition the white dwarf undergoes a simmer 
phase characterized by dipolar convection, with material flowing out of the center in a 
directed plume then circulating back in through the opposite side. The hottest points in 
the flow, and therefore where ignition occurs, are found on one side of the star displaced 
slightly from the center. As the simulations have been computed at much lower 
Reynolds number (-1,000) than the true star (~10 14 ), the real flow may be even less 
ordered, retaining a dipole flavor but contaminated by higher multipoles, making ignition 
a much more chaotic process. The geometry of ignition may also be influenced by the 
rotation rate of the star. In our models we therefore varied the number of ignition points, 
their distance from the center of the star, and the solid angle in which they were 
distributed. 

The physics of a putative transition to detonation, though still uncertain, has also 
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been elucidated in recent studies ' ' . As the deflagration flame propagates into lower 
density material, it thickens and slows to the point that turbulent eddies can mix hot fuel 
and cold ash, causing the burning rate to become highly irregular and potentially 
explosive. Like ignition, detonation should be a stochastic process which occurs only for 
the most extreme and intermittent turbulent energies, and may occur many times in 
different places. We therefore varied the criteria (critical Karlovitz number) for instituting 
a detonation (see supplementary material). 

Fig. 1 shows the final chemical structure (100 seconds after ignition) of the 
stellar debris for a sample of models. The inner regions of material, which were burned 
primarily in the turbulent, Rayleigh-Taylor unstable deflagration phase, consist of a 
patchy mixture of 56 Ni, intermediate mass elements (IMEs), and stable iron-group 
isotopes produced by electron capture (mostly 54 Fe and 58 Ni),. The subsequent detonation 
produced a smoother distribution of IMEs in the outer layers of the star. The synthetic 



model light curves, colors, and spectral time series agree very well with those of real 
events observed over the two months after explosion 24 (Figure 2), offering a strong 
validation of the model's predicted velocity structure and chemical stratification. On the 
other hand, it has been suggested that the mixing of electron capture elements throughout 
the inner layers - a feature generic to all multi-dimensional deflagration calculations - 
may be inconsistent with spectral observations of some SNe la taken at late phases . 

The models predict a range of 56 Ni yields, from 0.3 to 1.1 solar masses, depending 
on the initial conditions. Perhaps counter-intuitively, models with more robust ignition 
generally synthesize less 56 Ni. A strong ignition (i.e., numerous, symmetrically 
distributed sparks) increases the amount of burning in the deflagration phase during 
which the star expands and the density declines. This reduces the 56 Ni produced in the 
subsequent detonation wave. A weak or asymmetrical ignition, on the other hand, gives 
little deflagration burning and minimal pre-expansion, so that the detonation synthesizes 
abundant 56 Ni. For similar reasons, the 56 Ni mass also depended on the detonation 
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criteria, as noted in previous 1-D calculations . 

Given the range in 56 Ni masses, the model peak luminosities vary by a factor of 3, 
from 0.7 to 2.1 x 10 43 erg s" 1 which spans the range of normal SNe la, though does not 
reproduce the more extreme and peculiar sub- and super-luminous events . A 
substantial amount of the dispersion in SN la luminosities may therefore be the result of 
the stochastic and asymmetrical nature of the explosion itself, aside from any variations 
in the properties of the progenitor star. The rms dispersion in brightness depends on the 
wavelength band considered, being 45% in the blue, 27% in the red, and only 21% in the 
near infrared (-1.25 um). This reflects a self-regulating property of the radiative transfer 
— dimmer models are relatively cooler, and so radiate a greater percentage of their flux at 



longer wavelengths. The models thus confirm observational indications that SNe la are 
nearly standard (as opposed to merely standardizable) candles in the near infrared . 

Without introducing any artificial tuning, the luminosity of the models correlates 
with the light curve decline rate, giving a width-luminosity-relation (WLR) similar to that 
observed (Fig. 3). A correlation is also found between brightness and the color measured 
at peak. The rms scatter in the model WLR alone is 24%, while using both decline rate 
and color reduces the dispersion to 21% — similar to, but slightly greater than -16% that 

no 

is observed . The larger diversity seen in the models suggests that additional important 
physics may constrain the ignition and detonation conditions to a range narrower than 
that considered here. As in observations, the calibration of the models can be improved 
by using additional information from the light curve. For example, including the shape of 
the light curve in several optical and near infrared wavelength bands reduces the scatter 
to only 15%. 

For a given mass of 56 Ni, the residual scatter in the model WLR reflects 
individuating features of the supernova debris structure. The turbulent deflagration 
phase imprints density and chemical inhomogeneities which lead to variations in the time 
scale for photons to diffuse out of the debris. In addition, the global asymmetry - due to 
asymmetric ignition conditions or off-center detonation points - gives rise to anisotropic 
emission, so that the brightness and duration of most models vary by 20-30% depending 
of the angle from which they are viewed. While the adoption of a 2-D geometry may 
exaggerate global asymmetries, spectropolarization observations reveal that SNe la 
typically possess asphericity near or just below the level predicted here 14 . Dimmer 
supernovae tend to be more polarized, an observation consistent with our finding that 
dimmer models are more asymmetric due to relatively more burning in the turbulent 
deflagration phase. 



There are both theoretical suggestions and observational indications that the 
metallicity of the white dwarf will affect the 56 Ni yield at the -10% level. This is 
because the extra neutrons in trace elements such as 20 Ne lead to an increased synthesis 
of neutronized iron group elements ( 54 Fe and 58 Ni) at the expense of 56 Ni. To test the 
first-order effect of metallicity on the light curves, we varied the 56 Ni and metal 
abundances in the ejecta models according to predicted nucleosynthetic results. The 
resulting light curves (Fig. 4) show that both the peak luminosity and light curve duration 
decline with metallicity, in a manner roughly consistent with the WLR. Application of 
the WLR should therefore partially correct for metallicity variations, but with a residual 
error due to the different slope and normalization of the WLR at different metallicity. 
For extreme metallicity variations (from 0.3 to 3.0 times the solar value) the error can be 
as large as 4% in distance. The actual metallicity evolution over the range probed by 
cosmology experiments is much smaller than this, and we estimate that systematic errors 
in distance will likely be less than 2%. 

The models suggest that a substantial amount of the scatter in the observed WLR 
arises from the random sub-structures and viewing-angle effects that are predicted by 
multi-D explosions. In cosmological standard candle applications, these translate to easily 
reducible statistical errors. On the other hand, additional diversity arises from variations 
in metallicity and other properties of the progenitor star (e.g., carbon/oxygen ratio, central 
density, rotation) which may introduce a source of systematic error. Simulation offers one 
means to test how such variations influence the supernova light curve; for future dark 
energy experiments this will help to anticipate and limit errors arising from a shift in 
progenitor demographics over cosmic time. 
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Figure 1: Chemical structure of the ejected debris for a subset of the 
explosion models. Blue represents intermediate mass elements (i.e., sili- 
con, sulfur, calcium), green stable iron group elements produced by electron 
capture, and red 56 Ni. The turbulent inner regions reflect Rayleigh- Taylor 
and other instabilities that develop during the initial deflagration phase of 
burning. The subsequent detonation wave enhances the 56 Ni production in 
the center by burning remaining pockets of fuel. The lower density outer 
layers of debris, processed only by the detonation, consist of smoothly dis- 
tributed intermediate mass elements. 
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Figure 2: Synthetic multi-color light curves and spectra of a rep- 
resentative explosion model compared to observations of a nor- 
mal Type la supernova, a. The angle-averaged light curves of model 
DFD Jso_06_dc2 (solid lines) show good agreement with filtered observations 
of SN 20003du (Stanishev et a!., 2007; filled circles) in wavelength bands 
corresponding to the ultraviolet (U) optical (B,V,R), and near infrared (I), 
b. The synthetic spectra of the model (black lines) compare well to obser- 
vations of SN2003du (red lines) taken at times marked in days relative to 
B light curve maximum. Over time, as the remnant expands and thins, the 
spectral absorption features reflect the chemical composition of progressively 
deeper layers of debris, providing a strong test of the predicted compositional 
stratification of the model. 
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Figure 3: Correlation of the peak brightness of the models with 
their light curve duration and color. The sample includes 44 models 
each plotted for 30 different viewing angles. Solid circles denote models 
computed with the most likely range of detonation criteria, while open circles 
denote more extreme values, a. Relation between the peak brightness M B 
(measured in the logarithmic magnitude scale) and the light curve decline 
rate parameter AM 15 , defined as the decrease in B-band brightness from peak 
to 15 days after peak. The shaded band shows the approximate slope and 
spread of the observed width-luminosity relation, b. Relation between M B 
and the color parameter B-V measured at peak. The solid line shows the slope 
of the observed relation of Philips et al. (1999) but with the normalization 
shifted, as the models are systematically redder than observed SNe la by 7%, 
likely due to the approximate treatment of non-LTE effects. In observational 
studies, these two relations are usually fitted jointly as: M B = M Bfi + a(s — 
1) + (3[{B — ^)Bmax + (B — V)o], where s is a stretch parameter and (B — V)o 
is the color of a fiducial supernova. We take (B — V)o = and determine 
stretch using the first order relation: s — 1 — (AM15 — 1 . 1 ) / 1 . 7) . We find 
for the models fitted values of a = 2.25,/? = 4.45 and M b> q = —19.27 which 
are in reasonable agreement with those derived from the recent observational 
sample of Astier et al. (2006): a = 1.52,(3 = 1.57, and M Bfi = -19.31 + 
5 log 10 (f/o/70), where H is the Hubble parameter. 
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Figure 4: Effect of the metal content of the progenitor star popula- 
tion on the width-luminosity relation. The models explore two extreme 
values of the metallicity: 3 times (red points) and 0.3 times the solar value 
(blue points). For clarity, each model has been averaged over all viewing an- 
gles, and black lines connect similar explosion models of differing metallicity. 
The colored lines are linear fits to the width-luminosity relation of of the 
two metallicity samples separately. The diversity introduced by metallicity 
variations follows the general width-luminosity trend, but the slightly differ- 
ent normalization and slope of the relation for different metallicity samples 
indicates a potential source of systematic error in distance determinations. 
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Supplementary Material 



1. Explosion Simulations 
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For simulating the explosion process, we employed a code ' ' that treated the 
hydrodynamics in a higher-order Godunov scheme and followed the propagation of the 
thermonuclear flames in a level-set approach 34 . The effects of unresolved turbulence 
driving the propagation of the deflagration flame were accounted for by a subgrid-scale 
turbulence model 35 . This approach allowed for a self-consistent treatment of the flame 
propagation while avoiding tunable parameters. 

The deflagration level set was initiated representing a collection of spherical 
ignition sparks of radius 6 km distributed around the WD center, within a certain solid 
angle. From 20 to 150 ignition points were considered, randomly distributed from the 
center out to -300 km and in a solid angle whose opening varied from 60° to 360°. 

(see Table 1 for the ignition configurations of each model). Consequently, the ignition 
region forms a cone with the apex at the center of the WD. Inside this region, the sparks 
were randomly placed in angular direction. Their distribution in radial direction was 
randomly drawn from a Gaussian distribution. The apex angles of the cones and the 
standard deviations of the Gaussian distributions are listed in Table 1 . 

From these ignition configurations, the deflagration flame fronts evolved subject to 
buoyancy instabilities on large scales and driven by turbulent motions generated on 
smaller scales. The deflagration flames propagated from the center of the WD outwards 



and released energy, expanding the star. Until getting close to the surface of the WD, the 
deflagration flame is deformed on large scales by turbulence and thus its propagation is 
accelerated, but its microscopic structure remains "laminar", with its width and speed 
determined by radiative diffusion and conduction. Therefore it is stable against 
detonation. However, as the star expands and the flame burns towards its edge, the 
density directly ahead of the flame declines. The laminar flame structure becomes thicker 
and eventually small eddies can penetrate into the burning region and mix hot ash with 
cold fuel without immediately burning. This mixing process first begins when the 
Karlovitz number, Ka becomes greater than unity 36 , and it has long been speculated that a 
transition to detonation, were it to happen, would happen here . However, more recent 
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studies ' suggest that the first structures to form at the transition to distributed burning 
are too small to detonate. Greater mixing, and therefore a lower density and higher 
Karlovitz number are required. The largest mixed structures are formed when the eddy 
turnover time on the integral scale for the turbulence is equal to the nuclear burning time 
for the mixture, i.e., Damkohler number ~ 1. The effects of intermittency will raise this 
value some to perhaps Da ~ 10. Estimating the nuclear time scale with a small reaction 
network and assuming turbulent energies appropriate to the supernova 39 , this corresponds 
to Ka ~ 500. 

This possibility of a deflagration-to-detonation transition (DDT) was implemented in the 
models by determining Ka in each zone near the flame surface based upon the turbulent 
energy as derived from the subgrid-scale turbulence model. Here Ka is defined as Ka = 

1 10 

(6/Loib) where 6 is the laminar flame width, and the Gibson length, L G ib = 



[Siam/u'(Ax)] Ax, with Siam, the laminar flame speed and u', the turbulent rms fluctuation 
speed on the scale of the computational grid Ax. For a DDT to occur in our simulation, 
we required a minimal Karlovitz number and a certain range of fuel densities ahead of the 
flame, as given in Table 2. We considered five different values of critical Karlovitz 
number. Most models used a value in the most physically plausible range (Ka = 250, 
750) while other models explored more extreme values (Ka = 1, 1500, 2250). 

Once the flame fronts reached the critical conditions for detonation, a level set 
representing a detonation was initiated at the corresponding location. If other features of 
the flame reached critical DDT conditions, further detonations were initialized there on 
the basis of the same level set. Thus, a sequence of DDTs was possible and indeed 
realized in all models. The detonation level set was then propagated with an appropriate 
detonation velocity 40 , exhausting most of the remaining fuel. The level-set treatment of 
the detonation allowed us to prevent it from unphysically crossing ash regions left behind 
by the preceding deflagration phase. 

Each thermonuclear supernova explosion model is now specified by one of the 
deflagration ignitions configurations listed in Table 1 and one of the DDT conditions of 
Table 2. In total, 44 models were computed. 

Depending on the fuel density ahead of the flame, the C+O material of the WD was 
converted to a nuclear statistical equilibrium composition (NSE; modeled as a 
temperature-and density-dependent mixture of 56 Ni and alpha-particles), or to 



intermediate mass elements (modeled by a representative nucleus of A=30 and a nuclear 
binding energy of 8.17906 x 10 18 erg g" 1 ), or to oxygen. Electron captures in the NSE 
were accounted for in a parameterized way allowing us to differentiate radioactive 56 Ni 
from stable iron group elements. 

For select models, we carried out detailed post-processing nucleosynthesis calculations 
using Lagrangian tracer particles and a full nuclear network (see Figure 1 for two models 
discussed in the main text). For one of the more representative models 
(D2D_iso_06_dcl) we tabulated the detailed compositions comprising each of the coarse 
groups followed in the explosion calculations. These tables were used for the remaining 
models to interpolate the abundances of all elements, a prerequisite for the spectrum 
synthesis modeling described below. The light curves are less sensitive to the details of 
the abundance interpolation. 

The explosion simulations were run in cylindrical (r-z) geometry and a spatial resolution 
of 512x1024. Imposing rotational symmetry along the z-axis, the full star was represented 
in this two-dimensional setup. The hydrodynamical evolution of each model was 
followed for 100 seconds, well after burning had ceased and at which point the 
gravitational and internal energy densities were small (< 1%) relative to the kinetic 
energy density. At this time, the velocity structure was homologous (velocity 
proportional to radius) to better than a percent, indicating that the remnant had reached 
the phase of free-expansion. 



The explosion models all assumed the progenitor star had solar metallicity. To explore 
how varations in metallicity may influence the light curve, we changed the composition 
of the debris structure in each explosion model to reflect certain nucleosynthetic results. 
Timmes 29 showed in an analytical calculation that a higher metallicity leads to a greater 
production of stable iron group elements ( 54 Fe and 58 Ni) at the expense of 56 Ni. We 
therefore varied the 56 Ni mass according to M ni (Z) = M ni (Z=l) * (1 - 0.057 (Z-l)) , where 
the metallicity was changed to Z = 3 and Z = 0.3 times solar. The trace abundance of 
metals in unburned or partially burned material, which has some influence on the ejecta 
opacity, was also changed to reflect the progenitor metalicity. This approach does not 
account for the effect metallicity may have on the structure of the white dwarf, or on the 
dynamics of the explosion. 

2. Radiative Transfer Calculations 

Light curves and spectra of the models were computed using the multi- wavelength time- 
dependent radiative transfer code SEDONA 19 , which uses a Monte Carlo approach to 
solve the transport equation in arbitrary geometries. The final structure of each explosion 
model described above was remapped to a lower resolution (64x128) regular cylindrical 
grid to provide the initial conditions for the transfer code. SEDONA assumes the 
subsequent dynamics are given by homologous expansion and self-consistently calculates 
the temperature structure evolution by balancing sources of radiative heating and cooling. 

The light curves of SNe la are powered by the radioactive chain 56 Ni -> 56 Co -> 56 Fe. 
The decay releases primarily ~1 MeV gamma-rays, which deposit their energy in the 



ejecta mainly through Compton scattering and photo-electric absorption. The SEDONA 
code includes a detailed multi-wavelength transport scheme treating the emission, 
propagation, and absorption of gamma-rays. This provided the instantaneous rate and 
geometry of radioactive energy deposition, as well as predictions for the emergent 
gamma-ray light curves and spectra. 

Absorbed radioactive energy was assumed to be locally and instantaneously reprocessed 
into optical/UV photons, whose propagation was followed using a Monte Carlo method. 
Detailed non-grey opacities were applied, including the aggregate effects of over 10 
million bound-bound line transitions 41 treated in the expansion opacity formalism 42 . 
Atomic level populations were calculated assuming local thermodynamic equilibrium 
(LTE), typically a reasonable approximation for SNe la in the earlier epochs 43 . The 
radiation field, on the other hand, was not required to be in LTE, and an equivalent two- 
level equivalent atom (ETLA) formalism was used for the line source function: S = (1 - e) 
J + e B, where 8 is the ratio of absorptive opacity to total (scattering plus absorptive) 
opacity. A constant value 8 = 0.3 was used for all lines to approximate non-LTE effects 
based on comparison to previous line branching calculations. 

The most significant uncertainties in the light curve calculations relate to uncertainties in 
the calculation of the complex opacities/emissivities, in particular the likely inaccuracy 
and/or incompleteness of the atomic line database, and the limitations of the expansion 
opacity and ETLA formalisms. In the future, self-consistent multi-dimensional non-LTE 
calculations will help refine the quantitative accuracy of the models. 
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Table 1: Deflagration ignition configurations. The half apex angle of the cone centered 
on the WD in which the ignition is assumed to take place is denoted by a. 



Ignition setup 


Number of 


Minimal 


cos(a) 


Standard 




ignition 


distance of 




deviation 




kernels 


kernels 
[kernel radii] 




[km] 


DD2D_iso_01 


20 


1.0 


-1.0 


150.0 


DD20_iso_02 


50 


0.8 


-1.0 


150.0 


DD2D_iso_03 


60 


0.7 


-1.0 


150.0 


DD2D_iso_04 


80 


0.8 


-1.0 


150.0 


DD2D_iso_05 


90 


0.7 


-1.0 


150.0 


DD2D_iso_06 


100 


0.1 


-1.0 


150.0 


DD2D_iso_07 


100 


0.5 


-1.0 


150.0 


DD2D_iso_08 


150 


0.3 


-1.0 


150.0 


DD2D_asym_01 


120 


0.3 


-1.0 


75.0 


DD2D_asym_02 


105 


0.3 


-0.75 


75.0 


DD2D_asym_03 


90 


0.3 


-0.50 


75.0 


DD2D_asym_04 


75 


0.3 


-0.25 


75.0 


DD2D_asym_05 


60 


0.3 


0.0 


75.0 


DD2D_asym_06 


45 


0.3 


0.25 


75.0 


DD2D_asym_07 


30 


0.3 


0.50 


75.0 


DD2D_asym_08 


15 


0.3 


0.75 


75.0 



Table 2: Criteria for deflagration-to-detonation transitions 



DDT criterion (dc) 


Kctmin 


Pmin [10 7 g crrT 3 ] 


Pmax [10' g Cm"" 5 ] 


1 


1.0 


0.6 


1.75 


2 


250.0 


0.6 


1.20 


3 


750.0 


0.6 


1.20 


4 


1500.0 


0.6 


1.20 


5 


2250.0 


0.6 


1.20 



